<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN"
            "http://www.w3.org/TR/REC-html40/loose.dtd">
<HTML>
<HEAD>

<META http-equiv="Content-Type" content="text/html; charset=UTF-8">
<META name="GENERATOR" content="hevea 1.10">

<META name="Author" content="Julien Mairal">
<link rel="stylesheet" href="doc_spams.css">
<LINK rel="stylesheet" type="text/css" href="doc_spams.css">
<TITLE>Dictionary Learning and Matrix Factorization Toolbox</TITLE>
</HEAD>
<BODY >
<A HREF="doc_spams003.html"><IMG SRC="previous_motif.gif" ALT="Previous"></A>
<A HREF="index.html"><IMG SRC="contents_motif.gif" ALT="Up"></A>
<A HREF="doc_spams005.html"><IMG SRC="next_motif.gif" ALT="Next"></A>
<HR>
<H2 CLASS="section"><A NAME="htoc3">3</A>  Dictionary Learning and Matrix Factorization Toolbox</H2><UL>
<LI><A HREF="doc_spams004.html#toc1">Function mexTrainDL</A>
</LI><LI><A HREF="doc_spams004.html#toc2">Function mexTrainDL_Memory</A>
</LI><LI><A HREF="doc_spams004.html#toc3">Function nmf</A>
</LI><LI><A HREF="doc_spams004.html#toc4">Function nnsc</A>
</LI></UL>
<P>
This is the section for dictionary learning and matrix factorization, corresponding to [<A HREF="doc_spams009.html#mairal7">19</A>, <A HREF="doc_spams009.html#mairal9">20</A>].</P><H3 CLASS="subsection"><A NAME="toc1"></A><A NAME="htoc4">3.1</A>  Function mexTrainDL</H3><P>
This is the main function of the toolbox, implementing the learning algorithms of [<A HREF="doc_spams009.html#mairal9">20</A>]. 
Given a training set <I><B>x</B></I><SUP>1</SUP>,…, . It aims at solving
</P><TABLE CLASS="display dcenter"><TR VALIGN="middle"><TD CLASS="dcell">
</TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center">&nbsp;</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center">min</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center"><I><B>D</B></I> ∈ <FONT COLOR=red><I>C</I></FONT></TD></TR>
</TABLE></TD><TD CLASS="dcell"> </TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center">&nbsp;</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center">lim</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center"><I>n</I> → +∞</TD></TR>
</TABLE></TD><TD CLASS="dcell"> </TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center">1</TD></TR>
<TR><TD CLASS="hbar"></TD></TR>
<TR><TD CLASS="dcell" ALIGN="center"><I>n</I></TD></TR>
</TABLE></TD><TD CLASS="dcell"> </TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center"><I>n</I></TD></TR>
<TR><TD CLASS="dcell" ALIGN="center"><FONT SIZE=6>∑</FONT></TD></TR>
<TR><TD CLASS="dcell" ALIGN="center"><I>i</I>=1</TD></TR>
</TABLE></TD><TD CLASS="dcell"> </TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center">&nbsp;</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center">min</TD></TR>
<TR><TD CLASS="dcell" ALIGN="center">α<I><SUP>i</SUP></I></TD></TR>
</TABLE></TD><TD CLASS="dcell"> </TD><TD CLASS="dcell">⎛<BR>
⎜<BR>
⎝</TD><TD CLASS="dcell"><TABLE CLASS="display"><TR><TD CLASS="dcell" ALIGN="center">1</TD></TR>
<TR><TD CLASS="hbar"></TD></TR>
<TR><TD CLASS="dcell" ALIGN="center">2</TD></TR>
</TABLE></TD><TD CLASS="dcell"> ||<I><B>x</B><SUP>i</SUP></I>−<I><B>D</B></I>α<I><SUP>i</SUP></I>||<SUB>2</SUB><SUP>2</SUP> + ψ(α<I><SUP>i</SUP></I>)</TD><TD CLASS="dcell">⎞<BR>
⎟<BR>
⎠</TD><TD CLASS="dcell">.
    (1)</TD></TR>
</TABLE><P>
ψ is a sparsity-inducing regularizer and <FONT COLOR=red><I>C</I></FONT> is a constraint set for the dictionary. As shown in [<A HREF="doc_spams009.html#mairal9">20</A>] 
and in the help file below, various combinations can be used for ψ and <FONT COLOR=red><I>C</I></FONT> for solving different matrix factorization problems.
What is more, positivity constraints can be added to α as well. The function admits several modes for choosing the optimization parameters, using the parameter-free strategy proposed in [<A HREF="doc_spams009.html#mairal7">19</A>], or using the parameters <I>t</I><SUB>0</SUB> and ρ presented
in [<A HREF="doc_spams009.html#mairal9">20</A>]. <B>Note that for problems of a reasonable size, and when ψ is the ℓ<SUB>1</SUB>-norm, 
the function mexTrainDL_Memory can be faster but uses more memory.</B> </P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="#007F00">% <BR>
% Usage:   [D [model]]=mexTrainDL(X,param[,model]);<BR>
%          model is optional<BR>
%<BR>
% Name: mexTrainDL<BR>
%<BR>
% Description: mexTrainDL is an efficient implementation of the<BR>
%     dictionary learning technique presented in<BR>
%<BR>
%     "Online Learning for Matrix Factorization and Sparse Coding"<BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     arXiv:0908.0050<BR>
%     <BR>
%     "Online Dictionary Learning for Sparse Coding"      <BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     ICML 2009.<BR>
%<BR>
%     Note that if you use param.mode=1 or 2, if the training set has a<BR>
%     reasonable size and you have enough memory on your computer, you <BR>
%     should use mexTrainDL_Memory instead.<BR>
% <BR>
%<BR>
%     It addresses the dictionary learning problems<BR>
%        1) if param.mode=0<BR>
%     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2  s.t. ...<BR>
%                                                  ||alpha_i||_1 &lt;= lambda<BR>
%        2) if param.mode=1<BR>
%     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_1  s.t.  ...<BR>
%                                           ||x_i-Dalpha_i||_2^2 &lt;= lambda<BR>
%        3) if param.mode=2<BR>
%     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2 + ... <BR>
%                                  lambda||alpha_i||_1 + lambda_2||alpha_i||_2^2<BR>
%        4) if param.mode=3, the sparse coding is done with OMP<BR>
%     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2  s.t. ... <BR>
%                                                  ||alpha_i||_0 &lt;= lambda<BR>
%        5) if param.mode=4, the sparse coding is done with OMP<BR>
%     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_0  s.t.  ...<BR>
%                                           ||x_i-Dalpha_i||_2^2 &lt;= lambda<BR>
%        6) if param.mode=5, the sparse coding is done with OMP<BR>
%     min_{D in C} (1/n) sum_{i=1}^n 0.5||x_i-Dalpha_i||_2^2 +lambda||alpha_i||_0  <BR>
%                                           <BR>
%<BR>
%%     C is a convex set verifying<BR>
%        1) if param.modeD=0<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 &lt;= 1 }<BR>
%        2) if param.modeD=1<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... <BR>
%                                                  gamma1||d_j||_1 &lt;= 1 }<BR>
%        3) if param.modeD=2<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... <BR>
%                                  gamma1||d_j||_1 + gamma2 FL(d_j) &lt;= 1 }<BR>
%        4) if param.modeD=3<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  (1-gamma1)||d_j||_2^2 + ... <BR>
%                                  gamma1||d_j||_1 &lt;= 1 }<BR>
<BR>
%<BR>
%     Potentially, n can be very large with this algorithm.<BR>
%<BR>
% Inputs: X:  double m x n matrix   (input signals)<BR>
%               m is the signal size<BR>
%               n is the number of signals to decompose<BR>
%         param: struct<BR>
%            param.D: (optional) double m x p matrix   (dictionary)<BR>
%              p is the number of elements in the dictionary<BR>
%              When D is not provided, the dictionary is initialized <BR>
%              with random elements from the training set.<BR>
%           param.K (size of the dictionary, optional is param.D is provided)<BR>
%           param.lambda  (parameter)<BR>
%           param.lambda2  (optional, by default 0)<BR>
%           param.iter (number of iterations).  If a negative number is <BR>
%              provided it will perform the computation during the<BR>
%              corresponding number of seconds. For instance param.iter=-5<BR>
%              learns the dictionary during 5 seconds.<BR>
%           param.mode (optional, see above, by default 2) <BR>
%           param.posAlpha (optional, adds positivity constraints on the<BR>
%             coefficients, false by default, not compatible with <BR>
%             param.mode =3,4)<BR>
%           param.modeD (optional, see above, by default 0)<BR>
%           param.posD (optional, adds positivity constraints on the <BR>
%             dictionary, false by default, not compatible with <BR>
%             param.modeD=2)<BR>
%           param.gamma1 (optional parameter for param.modeD &gt;= 1)<BR>
%           param.gamma2 (optional parameter for param.modeD = 2)<BR>
%           param.batchsize (optional, size of the minibatch, by default <BR>
%              512)<BR>
%           param.iter_updateD (optional, number of BCD iterations for the dictionary<BR>
%              update step, by default 1)<BR>
%           param.modeParam (optimization mode).<BR>
%              1) if param.modeParam=0, the optimization uses the <BR>
%                 parameter free strategy of the ICML paper<BR>
%              2) if param.modeParam=1, the optimization uses the <BR>
%                 parameters rho as in arXiv:0908.0050<BR>
%              3) if param.modeParam=2, the optimization uses exponential <BR>
%                 decay weights with updates of the form <BR>
%                 A_{t} &lt;- rho A_{t-1} + alpha_t alpha_t^T<BR>
%           param.rho (optional) tuning parameter (see paper arXiv:0908.0050)<BR>
%           param.t0 (optional) tuning parameter (see paper arXiv:0908.0050)<BR>
%           param.clean (optional, true by default. prunes <BR>
%              automatically the dictionary from unused elements).<BR>
%           param.verbose (optional, true by default, increase verbosity)<BR>
%           param.numThreads (optional, number of threads for exploiting<BR>
%              multi-core / multi-cpus. By default, it takes the value -1,<BR>
%              which automatically selects all the available CPUs/cores).<BR>
%<BR>
% Output: <BR>
%         param.D: double m x p matrix   (dictionary)<BR>
%<BR>
% Note: this function admits a few experimental usages, which have not<BR>
%     been extensively tested:<BR>
%         - single precision setting <BR>
%<BR>
% Author: Julien Mairal, 2009</FONT></TD></TR>
</TABLE><P>The following piece of code contains usage examples:
</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="blue"><B>clear</B></FONT> <FONT COLOR="blue"><B>all</B></FONT>;<BR>
<BR>
I=double(imread(<FONT COLOR="red">'data/lena.png'</FONT>))/255;<BR>
<FONT COLOR="#007F00">% extract 8 x 8 patches</FONT><BR>
X=im2col(I,[8 8],<FONT COLOR="red">'sliding'</FONT>);<BR>
X=X-repmat(<FONT COLOR="blue"><B>mean</B></FONT>(X),[<FONT COLOR="blue"><B>size</B></FONT>(X,1) 1]);<BR>
X=X ./ repmat(<FONT COLOR="blue"><B>sqrt</B></FONT>(<FONT COLOR="blue"><B>sum</B></FONT>(X.^2)),[<FONT COLOR="blue"><B>size</B></FONT>(X,1) 1]);<BR>
<BR>
param.K=100;  <FONT COLOR="#007F00">% learns a dictionary with 100 elements</FONT><BR>
param.lambda=0.15;<BR>
param.numThreads=4; <FONT COLOR="#007F00">% number of threads</FONT><BR>
param.batchsize=400;<BR>
<BR>
param.iter=1000;  <FONT COLOR="#007F00">% let us see what happens after 1000 iterations.<BR>
<BR>
%%%%%%%%%% FIRST EXPERIMENT %%%%%%%%%%%</FONT><BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
D = mexTrainDL(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<BR>
param.approx=0;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
ImD=displayPatches(D);<BR>
<FONT COLOR="blue"><B>subplot</B></FONT>(1,3,1);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<FONT COLOR="blue"><B>drawnow</B></FONT>;<BR>
<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'*********** SECOND EXPERIMENT ***********\n'</FONT>);<BR>
<FONT COLOR="#007F00">%%%%%%%%%% SECOND EXPERIMENT %%%%%%%%%%%<BR>
% Train on half of the training set, then retrain on the second part</FONT><BR>
X1=X(:,1:<FONT COLOR="blue"><B>floor</B></FONT>(<FONT COLOR="blue"><B>size</B></FONT>(X,2)/2));<BR>
X2=X(:,<FONT COLOR="blue"><B>floor</B></FONT>(<FONT COLOR="blue"><B>size</B></FONT>(X,2)/2):<FONT COLOR="blue"><B>end</B></FONT>);<BR>
param.iter=500;<BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
[D model] = mexTrainDL(X1,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
<FONT COLOR="#007F00">% Then reuse the learned model to retrain a few iterations more.</FONT><BR>
param2=param;<BR>
param2.D=D;<BR>
[D model] = mexTrainDL(X2,param2,model);<BR>
<FONT COLOR="#007F00">%[D] = mexTrainDL(X,param);</FONT><BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<BR>
<FONT COLOR="#007F00">% let us add sparsity to the dictionary itself</FONT><BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'*********** THIRD EXPERIMENT ***********\n'</FONT>);<BR>
param.modeParam=0;<BR>
param.iter=1000;<BR>
param.gamma1=0.3;<BR>
param.modeD=1;<BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
[D] = mexTrainDL(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<FONT COLOR="blue"><B>tic<BR>
subplot</B></FONT>(1,3,2);<BR>
ImD=displayPatches(D);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>drawnow</B></FONT>;<BR>
<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'*********** FOURTH EXPERIMENT ***********\n'</FONT>);<BR>
param.modeParam=0;<BR>
param.iter=1000;<BR>
param.gamma1=0.3;<BR>
param.modeD=3;<BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
[D] = mexTrainDL(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<FONT COLOR="blue"><B>tic<BR>
subplot</B></FONT>(1,3,3);<BR>
ImD=displayPatches(D);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>drawnow</B></FONT>;</TD></TR>
</TABLE><H3 CLASS="subsection"><A NAME="toc2"></A><A NAME="htoc5">3.2</A>  Function mexTrainDL_Memory</H3><P>
Memory-consuming version of mexTrainDL. This function is well adapted to small/medium-size problems:
It requires storing all the coefficients α and is therefore impractical
for very large datasets. However, in many situations, one can afford this memory cost and it is better to use this method, which 
is faster than mexTrainDL.
Note that unlike mexTrainDL this function does not allow warm-restart.
</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="#007F00">% <BR>
% Usage:   [D]=mexTrainDL(X,param);<BR>
%<BR>
% Name: mexTrainDL_Memory<BR>
%<BR>
% Description: mexTrainDL_Memory is an efficient but memory consuming <BR>
%     variant of the dictionary learning technique presented in<BR>
%<BR>
%     "Online Learning for Matrix Factorization and Sparse Coding"<BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     arXiv:0908.0050<BR>
%     <BR>
%     "Online Dictionary Learning for Sparse Coding"      <BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     ICML 2009.<BR>
%<BR>
%     Contrary to the approaches above, the algorithm here <BR>
%        does require to store all the coefficients from all the training<BR>
%        signals. For this reason this variant can not be used with large<BR>
%        training sets, but is more efficient than the regular online<BR>
%        approach for training sets of reasonable size.<BR>
%<BR>
%     It addresses the dictionary learning problems<BR>
%        1) if param.mode=1<BR>
%     min_{D in C} (1/n) sum_{i=1}^n  ||alpha_i||_1  s.t.  ...<BR>
%                                         ||x_i-Dalpha_i||_2^2 &lt;= lambda<BR>
%        2) if param.mode=2<BR>
%     min_{D in C} (1/n) sum_{i=1}^n (1/2)||x_i-Dalpha_i||_2^2 + ... <BR>
%                                                      lambda||alpha_i||_1  <BR>
%<BR>
%     C is a convex set verifying<BR>
%        1) if param.modeD=0<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 &lt;= 1 }<BR>
%        1) if param.modeD=1<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... <BR>
%                                                  gamma1||d_j||_1 &lt;= 1 }<BR>
%        1) if param.modeD=2<BR>
%           C={  D in Real^{m x p}  s.t.  forall j,  ||d_j||_2^2 + ... <BR>
%                                  gamma1||d_j||_1 + gamma2 FL(d_j) &lt;= 1 }<BR>
%<BR>
%     Potentially, n can be very large with this algorithm.<BR>
%<BR>
% Inputs: X:  double m x n matrix   (input signals)<BR>
%               m is the signal size<BR>
%               n is the number of signals to decompose<BR>
%         param: struct<BR>
%            param.D: (optional) double m x p matrix   (dictionary)<BR>
%              p is the number of elements in the dictionary<BR>
%              When D is not provided, the dictionary is initialized <BR>
%              with random elements from the training set.<BR>
%           param.K (size of the dictionary, optional is param.D is provided)<BR>
%           param.lambda  (parameter)<BR>
%           param.iter (number of iterations).  If a negative number is <BR>
%              provided it will perform the computation during the<BR>
%              corresponding number of seconds. For instance param.iter=-5<BR>
%              learns the dictionary during 5 seconds.<BR>
%            param.mode (optional, see above, by default 2) <BR>
%            param.modeD (optional, see above, by default 0)<BR>
%            param.posD (optional, adds positivity constraints on the <BR>
%              dictionary, false by default, not compatible with <BR>
%              param.modeD=2)<BR>
%            param.gamma1 (optional parameter for param.modeD &gt;= 1)<BR>
%            param.gamma2 (optional parameter for param.modeD = 2)<BR>
%            param.batchsize (optional, size of the minibatch, by default <BR>
%              512)<BR>
%            param.iter_updateD (optional, number of BCD iterations for the dictionary <BR>
%                update step, by default 1)<BR>
%            param.modeParam (optimization mode).<BR>
%              1) if param.modeParam=0, the optimization uses the <BR>
%                 parameter free strategy of the ICML paper<BR>
%              2) if param.modeParam=1, the optimization uses the <BR>
%                 parameters rho as in arXiv:0908.0050<BR>
%              3) if param.modeParam=2, the optimization uses exponential <BR>
%                 decay weights with updates of the form <BR>
%                 A_{t} &lt;- rho A_{t-1} + alpha_t alpha_t^T<BR>
%            param.rho (optional) tuning parameter (see paper arXiv:0908.0050)<BR>
%            param.t0 (optional) tuning parameter (see paper arXiv:0908.0050)<BR>
%            param.clean (optional, true by default. prunes <BR>
%              automatically the dictionary from unused elements).<BR>
%            param.numThreads (optional, number of threads for exploiting<BR>
%              multi-core / multi-cpus. By default, it takes the value -1,<BR>
%              which automatically selects all the available CPUs/cores).<BR>
%<BR>
% Output: <BR>
%         param.D: double m x p matrix   (dictionary)<BR>
%<BR>
% Note: this function admits a few experimental usages, which have not<BR>
%     been extensively tested:<BR>
%         - single precision setting (even though the output alpha is double <BR>
%           precision)<BR>
%<BR>
% Author: Julien Mairal, 2009</FONT></TD></TR>
</TABLE><P>
The following piece of code contains usage examples:
</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="blue"><B>clear</B></FONT> <FONT COLOR="blue"><B>all</B></FONT>;<BR>
<BR>
I=double(imread(<FONT COLOR="red">'data/lena.png'</FONT>))/255;<BR>
<FONT COLOR="#007F00">% extract 8 x 8 patches</FONT><BR>
X=im2col(I,[8 8],<FONT COLOR="red">'sliding'</FONT>);<BR>
X=X-repmat(<FONT COLOR="blue"><B>mean</B></FONT>(X),[<FONT COLOR="blue"><B>size</B></FONT>(X,1) 1]);<BR>
X=X ./ repmat(<FONT COLOR="blue"><B>sqrt</B></FONT>(<FONT COLOR="blue"><B>sum</B></FONT>(X.^2)),[<FONT COLOR="blue"><B>size</B></FONT>(X,1) 1]);<BR>
X=X(:,1:10:<FONT COLOR="blue"><B>end</B></FONT>);<BR>
<BR>
param.K=200;  <FONT COLOR="#007F00">% learns a dictionary with 100 elements</FONT><BR>
param.lambda=0.15;<BR>
param.numThreads=4; <FONT COLOR="#007F00">% number of threads</FONT><BR>
<BR>
param.iter=100;  <FONT COLOR="#007F00">% let us see what happens after 100 iterations.<BR>
<BR>
%%%%%%%%%% FIRST EXPERIMENT %%%%%%%%%%%</FONT><BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
D = mexTrainDL_Memory(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
ImD=displayPatches(D);<BR>
<FONT COLOR="blue"><B>subplot</B></FONT>(1,3,1);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);<BR>
<BR>
<FONT COLOR="#007F00">%%%%%%%%%% SECOND EXPERIMENT %%%%%%%%%%%</FONT><BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
D = mexTrainDL(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
alpha=mexLasso(X,D,param);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-D*alpha).^2)+param.lambda*<FONT COLOR="blue"><B>sum</B></FONT>(<FONT COLOR="blue"><B>abs</B></FONT>(alpha)));<BR>
ImD=displayPatches(D);<BR>
<FONT COLOR="blue"><B>subplot</B></FONT>(1,3,2);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);</TD></TR>
</TABLE><H3 CLASS="subsection"><A NAME="toc3"></A><A NAME="htoc6">3.3</A>  Function nmf</H3><P>
This function is an example on how to use the function mexTrainDL for the
problem of non-negative matrix factorization formulated in [<A HREF="doc_spams009.html#lee2">16</A>]. Note
that mexTrainDL can be replaced by mexTrainDL_Memory in this function for
small or medium datasets.</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="#007F00">% <BR>
% Usage:   [U [,V]]=nmf(X,param);<BR>
%<BR>
% Name: nmf<BR>
%<BR>
% Description: mexTrainDL is an efficient implementation of the<BR>
%     non-negative matrix factorization technique presented in <BR>
%<BR>
%     "Online Learning for Matrix Factorization and Sparse Coding"<BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     arXiv:0908.0050<BR>
%     <BR>
%     "Online Dictionary Learning for Sparse Coding"      <BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     ICML 2009.<BR>
%<BR>
%     Potentially, n can be very large with this algorithm.<BR>
%<BR>
% Inputs: X:  double m x n matrix   (input signals)<BR>
%               m is the signal size<BR>
%               n is the number of signals to decompose<BR>
%         param: struct<BR>
%            param.K (number of required factors)<BR>
%            param.iter (number of iterations).  If a negative number <BR>
%              is provided it will perform the computation during the<BR>
%              corresponding number of seconds. For instance param.iter=-5<BR>
%              learns the dictionary during 5 seconds.<BR>
%            param.batchsize (optional, size of the minibatch, by default <BR>
%               512)<BR>
%            param.modeParam (optimization mode).<BR>
%               1) if param.modeParam=0, the optimization uses the <BR>
%                  parameter free strategy of the ICML paper<BR>
%               2) if param.modeParam=1, the optimization uses the <BR>
%                  parameters rho as in arXiv:0908.0050<BR>
%               3) if param.modeParam=2, the optimization uses exponential <BR>
%                  decay weights with updates of the form  <BR>
%                  A_{t} &lt;- rho A_{t-1} + alpha_t alpha_t^T<BR>
%            param.rho (optional) tuning parameter (see paper <BR>
%              arXiv:0908.0050)<BR>
%            param.t0 (optional) tuning parameter (see paper <BR>
%              arXiv:0908.0050)<BR>
%            param.clean (optional, true by default. prunes automatically <BR>
%              the dictionary from unused elements).<BR>
%            param.batch (optional, false by default, use batch learning <BR>
%              instead of online learning)<BR>
%            param.numThreads (optional, number of threads for exploiting<BR>
%                 multi-core / multi-cpus. By default, it takes the value -1,<BR>
%                 which automatically selects all the available CPUs/cores).<BR>
%         model: struct (optional) learned model for "retraining" the data.<BR>
%<BR>
% Output:<BR>
%         U: double m x p matrix   <BR>
%         V: double p x n matrix   (optional)<BR>
%         model: struct (optional) learned model to be used for <BR>
%           "retraining" the data.<BR>
%<BR>
% Author: Julien Mairal, 2009</FONT><BR>
<FONT COLOR="blue"><B>function</B></FONT> [U V] = nmf(X,param)<BR>
<BR>
param.lambda=0;<BR>
param.mode=2;<BR>
param.posAlpha=1;<BR>
param.posD=1;<BR>
param.whiten=0;<BR>
U=mexTrainDL(X,param);<BR>
param.pos=1;<BR>
<FONT COLOR="blue"><B>if</B></FONT> <FONT COLOR="blue"><B>nargout</B></FONT> == 2<BR>
   <FONT COLOR="blue"><B>if</B></FONT> <FONT COLOR="blue"><B>issparse</B></FONT>(X) <FONT COLOR="#007F00">% todo allow sparse matrices X for mexLasso</FONT><BR>
      maxbatch=<FONT COLOR="blue"><B>ceil</B></FONT>(10000000/<FONT COLOR="blue"><B>size</B></FONT>(X,1));<BR>
      <FONT COLOR="blue"><B>for</B></FONT> jj = 1:maxbatch:<FONT COLOR="blue"><B>size</B></FONT>(X,2)<BR>
         indbatch=jj:<FONT COLOR="blue"><B>min</B></FONT>((jj+maxbatch-1),<FONT COLOR="blue"><B>size</B></FONT>(X,2));<BR>
         Xb=<FONT COLOR="blue"><B>full</B></FONT>(X(:,indbatch));<BR>
         V(:,indbatch)=mexLasso(Xb,U,param);<BR>
      <FONT COLOR="blue"><B>end</B></FONT><BR>
   <FONT COLOR="blue"><B>else</B></FONT><BR>
      V=mexLasso(X,U,param);<BR>
   <FONT COLOR="blue"><B>end<BR>
end</B></FONT></TD></TR>
</TABLE><P>
The following piece of code contains usage examples:
</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="blue"><B>clear</B></FONT> <FONT COLOR="blue"><B>all</B></FONT>;<BR>
<BR>
I=double(imread(<FONT COLOR="red">'data/lena.png'</FONT>))/255;<BR>
<FONT COLOR="#007F00">% extract 8 x 8 patches</FONT><BR>
X=im2col(I,[16 16],<FONT COLOR="red">'sliding'</FONT>);<BR>
X=X(:,1:10:<FONT COLOR="blue"><B>end</B></FONT>);<BR>
X=X ./ repmat(<FONT COLOR="blue"><B>sqrt</B></FONT>(<FONT COLOR="blue"><B>sum</B></FONT>(X.^2)),[<FONT COLOR="blue"><B>size</B></FONT>(X,1) 1]);<BR>
<BR>
param.K=49;  <FONT COLOR="#007F00">% learns a dictionary with 100 elements</FONT><BR>
param.numThreads=4; <FONT COLOR="#007F00">% number of threads</FONT><BR>
<BR>
param.iter=-5;  <FONT COLOR="#007F00">% let us see what happens after 100 iterations.<BR>
<BR>
%%%%%%%%%% FIRST EXPERIMENT %%%%%%%%%%%</FONT><BR>
<FONT COLOR="blue"><B>tic</B></FONT><BR>
[U V] = nmf(X,param);<BR>
t=<FONT COLOR="blue"><B>toc</B></FONT>;<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'time of computation for Dictionary Learning: %f\n'</FONT>,t);<BR>
<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'Evaluating cost function...\n'</FONT>);<BR>
R=<FONT COLOR="blue"><B>mean</B></FONT>(0.5*<FONT COLOR="blue"><B>sum</B></FONT>((X-U*V).^2));<BR>
ImD=displayPatches(U);<BR>
<FONT COLOR="blue"><B>imagesc</B></FONT>(ImD); <FONT COLOR="blue"><B>colormap</B></FONT>(<FONT COLOR="red">'gray'</FONT>);<BR>
<FONT COLOR="blue"><B>fprintf</B></FONT>(<FONT COLOR="red">'objective function: %f\n'</FONT>,R);</TD></TR>
</TABLE><H3 CLASS="subsection"><A NAME="toc4"></A><A NAME="htoc7">3.4</A>  Function nnsc</H3><P>
This function is an example on how to use the function mexTrainDL for the
problem of non-negative sparse coding as defined in [<A HREF="doc_spams009.html#hoyer">13</A>]. Note that
mexTrainDL can be replaced by mexTrainDL_Memory in this function for small or
medium datasets.</P><TABLE CLASS="lstframe" STYLE="padding:1ex;background-color:#E5EDFF;"><TR><TD CLASS="mouselstlisting"><FONT COLOR="#007F00">% <BR>
% Usage:   [U [,V]]=nnsc(X,param);<BR>
%<BR>
% Name: nmf<BR>
%<BR>
% Description: mexTrainDL is an efficient implementation of the<BR>
%     non-negative sparse coding technique presented in <BR>
%<BR>
%     "Online Learning for Matrix Factorization and Sparse Coding"<BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     arXiv:0908.0050<BR>
%     <BR>
%     "Online Dictionary Learning for Sparse Coding"      <BR>
%     by Julien Mairal, Francis Bach, Jean Ponce and Guillermo Sapiro<BR>
%     ICML 2009.<BR>
%<BR>
%     Potentially, n can be very large with this algorithm.<BR>
%<BR>
% Inputs: X:  double m x n matrix   (input signals)<BR>
%               m is the signal size<BR>
%               n is the number of signals to decompose<BR>
%         param: struct<BR>
%           param.K (number of required factors)<BR>
%           param.lambda (parameter)<BR>
%           param.iter (number of iterations).  If a negative number <BR>
%              is provided it will perform the computation during the<BR>
%              corresponding number of seconds. For instance param.iter=-5<BR>
%              learns the dictionary during 5 seconds.<BR>
%           param.batchsize (optional, size of the minibatch, by default <BR>
%              512)<BR>
%           param.modeParam (optimization mode).<BR>
%              1) if param.modeParam=0, the optimization uses the <BR>
%                 parameter free strategy of the ICML paper<BR>
%              2) if param.modeParam=1, the optimization uses the <BR>
%                 parameters rho as in arXiv:0908.0050<BR>
%              3) if param.modeParam=2, the optimization uses exponential <BR>
%                 decay weights with updates of the form <BR>
%                 A_{t} &lt;- rho A_{t-1} + alpha_t alpha_t^T<BR>
%           param.rho (optional) tuning parameter (see paper<BR>
%              arXiv:0908.0050)<BR>
%           param.t0 (optional) tuning parameter (see paper <BR>
%              arXiv:0908.0050)<BR>
%           param.clean (optional, true by default. prunes automatically <BR>
%              the dictionary from unused elements).<BR>
%           param.batch (optional, false by default, use batch learning <BR>
%              instead of online learning)<BR>
%           param.numThreads (optional, number of threads for exploiting<BR>
%              multi-core / multi-cpus. By default, it takes the value -1,<BR>
%              which automatically selects all the available CPUs/cores).<BR>
%         model: struct (optional) learned model for "retraining" the data.<BR>
%<BR>
% Output:<BR>
%         U: double m x p matrix   <BR>
%         V: double p x n matrix   (optional)<BR>
%         model: struct (optional) learned model to be used for <BR>
%            "retraining" the data.<BR>
%<BR>
% Author: Julien Mairal, 2009</FONT><BR>
<BR>
[U V] = <FONT COLOR="blue"><B>function</B></FONT> nnsc(X,param)<BR>
<BR>
param.mode=2;<BR>
param.posAlpha=1;<BR>
param.posD=1;<BR>
param.whiten=0;<BR>
U=mexTrainDL(X,param);<BR>
param.pos=1;<BR>
<FONT COLOR="blue"><B>if</B></FONT> <FONT COLOR="blue"><B>nargout</B></FONT> == 2<BR>
   V=mexLasso(X,U,param);<BR>
<FONT COLOR="blue"><B>end</B></FONT></TD></TR>
</TABLE><HR>
<A HREF="doc_spams003.html"><IMG SRC="previous_motif.gif" ALT="Previous"></A>
<A HREF="index.html"><IMG SRC="contents_motif.gif" ALT="Up"></A>
<A HREF="doc_spams005.html"><IMG SRC="next_motif.gif" ALT="Next"></A>
</BODY>
</HTML>
